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ABSTRACT 

Recent observations show that supermassive black holes (SMBHs) with ~ 10 9 Mg 
exist at redshift z > 6; how they form has yet to be explained. A promising formation 
channel is the so-called direct collapse model, which posits that a massive seed BH 
forms through gravitational collapse of a ~ 10 5 M 0 supermassive star. We study the 
evolution of such a supermassive star growing by rapid mass accretion. The inter¬ 
nal stellar structure is also followed consistently in our calculation. In particular, we 
examine the impact of time-dependent mass accretion of repeating burst and quies¬ 
cent phases that are expected to occur with a self-gravitating circumstellar disk. We 
show that the stellar evolution with such episodic accretion differs qualitatively from 
that expected with a constant accretion rate, even if the mean accretion rate is the 
same. Unlike the case of constant mass accretion, whereby the star expands roughly 
following I?* ~ 2.6 x 10 3 .Rg(M*/100 Mg) 1 / 2 , the protostar can substantially contract 
during the quiescent phases between accretion bursts. The stellar effective temperature 
and ionizing photon emissivity increase accordingly as the star contracts, which can 
cause strong ionizing feedback and halt the mass accretion onto the star. With a fixed 
duration of the quiescent phase A t q , such contraction occurs in early evolutionary 
phases, i.e. for M* < 10 3 Mg with At q ~ 10 3 yr. For later epochs and larger masses 
but the same At q , contraction is negligible even during quiescent phases. With larger 
quiescent times At q , however, the star continues to contract during quiescent phases 
even for the higher stellar masses. We show that such behavior is well understood 
by comparing the interval time and the thermal relaxation time for a bloated surface 
layer. We conclude that the UV radiative feedback becomes effective if the quiescent 
phase associated by the burst accretion is longer than ~ 10 3 yr, which is possible in 
an accretion disk forming in the direct collapse model. 

Key words: cosmology: theory, early Universe, stars: formation, galaxies: formation, 
quasars: supermassive black holes 


1 INTRODUCTION 


Recent observations have revealed the existence of super¬ 
massive black holes (SMBHs) with masses of ~ 10 9 M 0 at 
x > 6 (e.g., |Wu et al.||2015| |Marziani fe Sulentic||2012| 
Mortlock et al. 201 1| ). The origin and the rapid growth of 
the early SMBHs remain to be elucidated. 

Possible seeds of the SMBHs are remnant BHs of Pop¬ 
ulation III stars with Mbh ~ 100 Mq (e.g., Madau & Rees 


2001 Schneider et al. 2002). Although such a ~ 100 M 0 


seed BH can just barely attain a mass of ~ 10 J M 0 by 
the epoch of z — 6, by accreting material at the Eddington 
rate, recent studies suggest difficulties in this model. For in¬ 


stance, radiative feedback from a BH accretion disk easily 
suppresses the gas supply from the intergalactic medium, so 
that the accretion rates fall far below the Eddington values 
(e.g., Alvarez, Wise & Abel 2009 Jeon et al. 20121. The 


BH growth time then becomes much longer than the age of 
the universe at z ~ 6, 


An alternative model which circumvents these difficul¬ 
ties is the so-called direct collapse model. The model as¬ 
sumes that SMBHs are built from larger seed BHs with 
Mbh ~ 10 5 Mq that are formed directly by gravitational 
collapse of supermassive stars (SMSs) of similar masses (e.g., 
|Bromm fe Loeb||2003| ). The BH growth time is significantly 
shortened in this case; a ~ 10 5 M 0 seed BH can easily grow 
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to become a ~ 10 J Mg SMBH by accretion at the Edding¬ 
ton rate. At reduced accretion rates, e.g. due to radiative 


lapse model. In this case, a protostar continues to expand 
and does not contract even after the stellar mass greatly ex- 


feedback from a BH-disk system, a ~ 10 9 Mg SMBH can 

ceeds 100 M© (Hosokawa, Omukai & Yorke|2012) 

The stel- 

still form at z ~ 6 if the mean accretion rate is higher than 

lar radius reaches ~ 10“* Rq for M* > lO 1 * Mg l 

Hosokawa 


about 50 % of MBdd- Indeed, cosmological simulations show 
that this is possible with efficient cold accretion flows re¬ 


alized in the epoch of the first-galaxy formation (e.g., Di 
Matteo et al.|[2012 |. 


A crucial assumption in the direct collapse model is that 
sufficiently massive stars are formed in a gas cloud. SMSs 
are supposed to form i n so-called atomi c-coo ling haloes with 
T vir 


10 4 K(e.g., Inayoshi fe Qmukai||2012 Agarwal et al. 


2014 Visbal, Haiman & Bryan 20141. If molecular hydro¬ 


gen cooling is inhibited in these massive primordial haloes 
by strong ultra-violet radiation, a gas cloud gravitationally 
collapses via atomic hydrogen cooling nearly isothermally 


at T ~ 8000 K (e.g., Omukai 2001). An embryo protostar 
eventually forms and begins to grow via gas accretion from a 


Van Borm et al.]2014i, analogous to normal Population III the star (e.g., Greif et al. 2012), the accretion rate dras 


surrounding envelope (e.g., Inayoshi, Omukai & Tasker 2014 


tostar” has a low effective temperature of Teff — 5000 K 
and the ionizing photon emissivity for M* < 10 4 Mg is only 
< 10 45 s , several orders of magnitude lower than for main- 
sequence SMSs. The resulting UV feedback should be too 
weak to halt the mass accretion. The rapid mass accretion 
should thus continue, so that SMSs finally form. 

The previous studies consider only constant accre¬ 
tion rates. In more realistic situations, however, the mass 
accretion onto a growing SMS should be dynamic with 
highly time-dependent accretion histories. For normal Pop 
III cases, for instance, an accretion disk formed around a 
protostar becomes gravitationally unstable and fragments 
(e.g., Stacy, Greif fe Bromm|2010 Greif et al.|2011 1. When 
such fragments migrate inward in the disk and accrete onto 


star formation (|Yoshida, Omukai fe Hernquist 20081. The tically increases causing so-called “burst accretion” (e.g., 


accretion rate at this stage has the well-known temperature 
dependence 


Vorobyov, DeSouza & Basu 2013). In this case, the accre- 
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tion history is roughly divided into two phases: the burst 
phases with M* ~ 0.1 Mg yr -1 and quiescent phases with 
M* ~ 10 -4 Mg yr -1 , whose durations are < 10 2 yr and 
~ 10 3 — 10 4 yr. Such episodic accretion is also expected to 


where c s is the sound speed and G is the gravitational con¬ 
stant. Because of the higher gas temperature during the col¬ 
lapse stage, the accretion rate is much higher than that for 
normal Pop III cases. If such rapid mass accretion continues 
for ~ 1 Myr, which is the lifetime of massive stars, SMSs 
with 10° — 10 6 Mg can form and ultimately provide massive 
BH seeds after their deaths. 

A possible obstacle for the growth of SMSs via rapid 
mass accretion is strong radiative feedback from the central 
massive protostar. Even when a normal Pop III star with 
~ 100 Mg forms by accretion, stellar UV radiation becomes 
so strong to create an Hu region that dynamically expands 


occur in the direct collapse model (e.g., Inayoshi & Haiman 


2014). Indeed, recent high-resolution numerical simulations 


report signatures of the disk fragmentation in atomic-cooling 


through the accretion envelope(e.g. McKee & Tan 20081. 


The gas accretion is then halted; the final stellar mass is 


effectively determined by this mechanism (Hosokawa et al. 


2011 Hirano et al. 20141. For the case of SMS formation 


considered here, similar, or even stronger UV feedback might 
significantly reduce or even halt the gas mass accretion. This 
could occur because a SMS should emit copious amounts of 
photons; the stellar luminosity is nearly at the Eddington 
value, which for massive stars is proportional to stellar mass. 
The resulting UV feedback could thus be much stronger than 
for normal Pop III cases. 

The strength of the stellar radiative feedback depends 
critically on the ionizing photon emissivity of primordial pro¬ 
tostars, which is determined by their evolution while accret¬ 
ing. For normal Pop III cases with mean accretion rates of 
a few xl0~ 3 Mg yr -1 , an accreting protostar has a large 
radius during the early evolution but eventually contracts 
to the zero-age main-sequence (ZAMS) for M» ~ 100 Mg 
(e.g., O mukai fe Palla|20 03|. When the protostar contracts, 
its ionizing photon emissivity increases and UV radiative 
feedback becomes effective. 

By contrast, the protostellar evolution should qualita¬ 
tively differ from the above at very rapid accretion rates 
M* > M*, cr = 4 x 10 -2 M 0 yr -1 expected in the direct col- 


haloes (Regan, Johansson & Haehnelt 2014 Becerra et al. 
[20l5] ). 

In this paper, we study the evolution of SMSs with time- 
dependent accretion histories in which a number of accre¬ 
tion bursts occur. Whereas the radiative feedback is gener¬ 
ally weak for constant accretion, a time-dependent accretion 
rate might significantly affect the protostellar evolution and 
hence the feedback strength. In quiescent phases between 
the bursts, for instance, accretion rates should fall below the 
critical rate M*, cr = 4 x 10~ 2 Mg yr -1 , above which a proto¬ 
star enters the supergiant stage. The protostar may contract 
in this case, so that the effective temperature rises and ra¬ 
diative feedback becomes effective. Here, we show that, by 
solving the stellar interior structure numerically, this should 
occur even with a mean accretion rate of ~ 0.1 Mg yr -1 , if 
the quiescent phase lasts for periods > 10 3 yr. 

The remainder of this paper is organized as follows. In 
Section [2] our numerical methods and modeling of the burst 
accretion are explained. The numerical results are presented 
in Section [3J where we also compare the protostellar evolu¬ 
tion for constant accretion rates to the evolution for time- 
dependent rates with a number of accretion bursts. We dis¬ 
cuss implications of these results in Section[4]and summarize 
our conclusions in Section [5] 


2 NUMERICAL METHOD AND MODELING 
OF BURST ACCRETION 

2.1 Numerical method 

To calculate the stellar evolution, we use the numerical code 


originally developed by Yorke & Bodenheimer (20081 with 
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additional improvements as in Hosokawa et al. (2013). The 


code solves the following four basic equations of the stellar 
evolution: 


dr 

dm 

dP 

dm 


4-7T r 2 p 

Gm 

4nr 4 


dl_ dT 8 dP 

dm nuc Cp dt p dt 


dT 

dm 


GmT 

4nr 4 P 


V, 


( 2 ) 

(3) 

(4) 

(5) 


where m is the mass coordinate, r the radial distance from 
the center, P the total (radiation plus gas) pressure, l the 
local luminosity, T the temperature, p the density, E nuc the 
net energy generation rate by nuclear fusion, cp the isobaric 
specific heat, V = <?lnT/91nP the actual temperature gra¬ 
dient, and 8 = — (dlnp/dlnT) p . The temperature gradi¬ 
ent V is evaluated in the context of mixing-length theory 
for convective layers. The nuclear reactions of hydrogen and 
helium burning are included in the code. 

In Eq. @, hydrostatic equilibrium is always assumed 
and the inertial term is omitted. This assumption is well- 
founded, because the local free-fall time in the stellar in¬ 
terior is generally much shorter than 10 years, the typical 
timescale over which the mass accretion rate varies (also see 


Sec. 2.2 below). Note that we consider the time-derivative 
in Eq. (|4j). Thus, thermal equilibrium is not assumed in our 
calculations. 

Mass accretion is implemented by adding a mass M* At 
at each time step to the outermost layer of stellar models. We 
assume that the newly accreting gas has the same physical 
quantities as in the stellar atmosphere. This will be realized 
if the accreting materials have enough time to adjust their 
thermal states to those in the atmosphere, e.g., when slowly 
approaching the stellar surface orbiting in a circumstellar 
disk. This is obviously a limiting case, and the accreting 
gas may have some additional thermal energy. We model 
this potential difference using a free-parameter r/, which is 
the fraction of the energy advected into the star with the 
accreting gas to the released gravitational energy, 


V = 


= L, 


/ GM, M, 


R, 


( 6 ) 


where L*, acc is the additional energy input to the star. The 
values of rj adopted in our set of models are summarized 
in Ta ble [l| ( also see Section 2.2 I. As explained in Hosokawa 
et al. (20131, however, rj actually has little effect on the stel¬ 


lar evolution except during early stages when M* <100 Mg. 

We follow the growth of a primordial supermassive pro¬ 
tostar by numerically solving the stellar structure equations. 
We consider several models of accreting protostars that grow 


via episodic accretion (see 2.2 below). Our calculations be¬ 


gin with an initial pre-main sequence stellar model of mass 
2 Mq and radius 25 Rq , whereby the initial interior struc¬ 
ture is determined by solving the Lane-Emden equation with 
a polytropic index n = 1.5. We assume that the initial model 
and newly accreting gas have the primordial composition 
with X = 0.72 and Y = 0.28. 


2.2 Modeling of burst accretion 

We model accretion histories expected in the case of a 
self-gravitating accretion disk. Numerical simulations show 
that, in such cases, the accretion rate sometimes sharply 
increases like a burst. For instance, such a burst occurs 
when a fragment formed via the gravitational instability mi¬ 
grates inward through the disk to be accreted onto the star 
I Vorobyov, DeSouza fe Basu||2013 ). Normally, such a burst 
is followed by a quiescent phase, when the mass accretion al¬ 
most ceases. Another accretion burst can be triggered when 
the disk self-gravity is effective again due to mass growth of 
the disk from a surrounding envelope. We study the effect 
of burst accretion for the case of the direct collapse model. 

We model accretion histories with periodic accretion 
bursts (e.g., Baraffe, Chabrier & Gallardo 2009 Hosokawa, 


Offner & K rumholz|2011 ), by assuming alternating constant 

high and low accretion rates for the burst and quiescent 
phases, respectively: b and M*, q . The durations of these 

phases are parametrized with Atb and At q . We also assume 
that the transition from a burst to a quiescent phase takes 
place over a finite duration Aft- The accretion rates dur¬ 
ing the transition phases are given by linear interpolation. 
In total, we have five parameters that characterize a model 
accretion history: b , M*, q , Atb, At q , and At t . 

We examine four different cases with different sets of 
these parameters as listed in Table [T] The values are cho¬ 
sen so that the mean accretion rate is 0.1 Mg/yr, which 
is expected for the direct collapse model. We compare the 
results of the four cases with effectively the same mean ac¬ 
cretion rate. We start the evolutionary calculations with an 
accretion burst. 

The finite transition phase is needed for numerical sta¬ 
bility when calculating stellar evolution with highly time- 
dependent mass accretion histories. We assume finite tran¬ 
sition phases with Atb < At t < At q (see Table [l] for the 
adopted values). We have checked that varying the transi¬ 
tion time by a factor of a few has little effect on the results. 


3 RESULTS 

3.1 Evolution with constant accretion rates 

3.1.1 Normal Pop III case with 10 -3 Mg yr -1 

Before presenting the protostellar evolution with the burst 
accretion, we briefly describe the fiducial case with a con¬ 
stant accretion rate for comparison. In the typical case 
with the accretion rate of 1O~ 3 M 0 yr _1 , a protostar goes 
through several distinct evolutionary stages. Let us consider 


the following two timescales (e.g., Stahler, Palla & Salpeter 


1986; Omukai & Pall a|2003 Hosokawa & Omukai 20091: the 
Kelvin-Helnrholtz (KH) time 


Ikh = 


GMl 

R,L, 


(7) 


which is the thermal relaxation time over which the grav¬ 
itational energy is released by radiation, and the accretion 
time 


tacc — 


M* 

M, 


( 8 ) 
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Table 1. Models of episodic accretion whose mean accretion rates are about 0.1 A/g./yr. 


Model 

A 

B 

C 

D 

Duration of the burst phase Atb [ yr ] 

25 

50 

100 

500 

Duration of the quiescent phase At q [ yr ] 

270 

540 

1080 

5400 

Accretion rate in the burst phase M * b [ Mq /yr ] 

1 

1 

1 

1 

Accretion rate in the quiescent phase M*, q [ M@/y r ] 

10“ 3 

10“ 3 

10“ 3 

10~ 3 

Transition time Att [ yr ] 

50 

100 

200 

1000 

Radiative efficiency 77 

0.1 

0.1 

0.1 

0.01 


which is the characteristic stellar growth time. In the earliest 
phase when M* < 5 Mg, the star evolves adiabatically, satis¬ 
fying t acc < tKH- As the stellar luminosity L* increases with 
increasing stellar mass, the KH time becomes shorter. The 
luminosity increases because the the stellar interior opacity 
decreases according to Kramers’ law (k oc pT~ i b ) when the 
temperature increases. As the opacity decreases, the heat 
accumulated in the interior gradually escapes outward by 
radiative heat transport. The timescale balance is finally in¬ 
verted and the protostar begins to contract radiating the 
energy away from the surface. The KH contraction stage 
begins when M* > 8 Mg in the case of M = 10 -3 Mq yr -1 , 
as can be seen in Figure [l] The stellar interior temperature 
rises during the KH contraction stage, and finally hydrogen 
burning begins at the center when M* ~ 40 Mq. After this 
point, the evolution of the stellar radius closely traces the 
mass-radius relation of a ZAMS star. 

Figure [I] shows that the rapid rise of the stellar ionizing 
photon emissivity is synchronous with the KH contraction. 
Before KH contraction begins, the ionizing photon emissiv¬ 
ity is only ~ 10 3 ' sec' 1 . As the star contracts and the ef¬ 
fective temperature rises, however, the emissivity quickly 
reaches > 10 49 sec -1 . Radiation hydrodynamic simulations 
show that the stellar radiation feedback operates in the late 
KH stage and eventually shuts off the mass accretion onto 
the star (e.g., Hosokawa et al.|2011 1. 


3.1.2 Direct collapse case with 0.1 Mq yr 1 

Here, we consider protostellar evolution at the high con¬ 
stant accretion rate of 0.1 Mq yr -1 . Unlike the case for 
10 -3 Mq yr -1 , the protostar does not evolve to the KH 
contraction stage (see the dashed lines in Fig. [lj. The star 
instead continues to expand nearly monotonically with in¬ 
creasing mass. 

The timescale inversion described in Section 13.1.11 oc¬ 
curs at M* ~ 22 Mq ( t ~ 200 yr) in this case (Fig. [2j|. 
The star evolves adiabatically before this epoch. Since the 
rapid accretion enhances the average entropy in the stellar 
interior, the stellar radius at a given mass is larger than for 
10 -3 Mq yr -1 . 

The protostar continues expanding even when Ikh < 
face; the star radiates the internal energy through the sur¬ 
face, which normally makes the star contract. (Hosokawa, | 
Omukai & Yorke (2012) show that most of the stellar inte¬ 


rior actually contracts even in this bloating phase. Basically, 
only the newly accreted surface layer inflates by absorbing 
the outward heat flux coming from the contracting interior. 


The mass of the bloating surface layer has a small fraction 
of the total stellar mass. In the surface layer, the opacity 
is mostly due to H - bound-free absorption that has a very 
strong temperature-dependence. As a result, the stellar ef¬ 
fective temperature is locked at a constant value ~ 5000 K 
as in the case of red-giants. From relations T e fr ~ 5000 K 
and L* = 47ri?*T 4 ff ~ ^Edd, we can easily show 


R t ~ 2.6 x 10 3 Rq 


M * 


100 Mq 


1/2 


(9) 


which explains the numerical results very well. The evolu¬ 
tion of ionizing photon emissivity also differs from the case 
with M = 10 -3 Mq yr -1 . Even when M* > 100 Mq, the 
protostar still has very low ionizing photon emissivity be¬ 
cause of the low effective temperature. The resulting stellar 
UV feedback is then too weak to halt the rapid mass accre¬ 
tion (see Section [4|. 

We have found that in spite of the timescale imbalance 
Ikh < tacc, the star continues to expand. The apparent dis¬ 
crepancy is explained by the fact that the normal definition 
of f K H (Eq. 0 is a global relationship and does not consider 
the stellar interior structure; as we discussed above, most 
of the stellar mass is concentrated near the center and the 
bloated surface layer contains only a small fraction of the 
stellar mass (also see Figure 2 in Hosokawa et al. (201 31. 

We therefore must evaluate the local KH timescale for 
the surface layer by considering the actual mass distribution 
within a star, 


Ikh, surf — 


/ / SradTdm 

7dl ’ 


( 10 ) 


where s ra d is the entropy of radiation and / is a dimen¬ 
sionless constant of 0(1). The latter factor is introduced to 
represent that only a fraction / of the total entropy is car¬ 
ried away over the timescale. We compare the above local 
KH time with the local accretion timescale for the surface 
layer, 


^acc.surf 


//dm 

M* 


( 11 ) 


For consistency we also include the same factor / in Eq. 
(111. In the following, we set / = 0.4 as a fit to our nu¬ 
merical results. We take the integration range in Eqs. (101 
and as 0.7 M» ^ m ^ M, to cover the surface layer. 
We have checked that varying the lower bound does not sig¬ 
nificantly affect the main results. Figure [2] shows that the 
local accretion time t acc ,surf still remains shorter than the 
corresponding KH time Ikh, surf even when the global KH 
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timescale satisfies i ac c > £kh- The star continues to expand 
because, with the rapid mass accretion, the gas accumulates 
in the surface layer faster than it can thermally relax. In 
other words, part of the outgoing heat flux is trapped by 
the accreting gas, and the inflated surface layer cannot con¬ 
tract. 

Below we show that these two local timescales f acCiSur f 
and fKH.surf are also useful for understanding stellar evolu¬ 
tion for the case of burst accretion. 


3.2 Evolution with burst accretion 

3.2.1 Model C 


We first focus on our Model C as a fiducial case for stellar 
evolution with burst accretion. In this case, the burst and 
quiescent phases have the accretion rates M*,b = 1 Mq yr -1 
and M*, q = 10 -3 Mq yr -1 for the durations of Aib = 100 
yr and At q = 1080 yr (Table 1). The protostellar evolu¬ 
tion is shown in Figure [3] where we note substantial dif¬ 
ferences in the evolution from the case with the constant 
rate 0.1 Mq yr -1 . First of all, the star contracts during the 
quiescent phases. At 300 < t (yr) < 1500, for instance, the 
stellar radius decreases to 100 Rq, which is more than 10 
times smaller than the super-giant protostar for the case 
M = 0.1 Mq yr -1 at the same mass. When the star con¬ 
tracts, the ionizing photon emissivity increases and becomes 
comparable to that of the corresponding ZAMS star. 

Because the accretion rate in the quiescent phase is 
lower than the critical rate for maintaining bloating of the 
star, ~ 4 x 10 -2 Mq yr -1 , the star contracts as in the nor¬ 
mal KH contraction stage. To see this, we analytically derive 
an equation which describes the time evolution of the radius 
assuming that the contraction occurs over the KH timescale, 

R* 1 c t -to 1 

Rq \R*, q /Rq 1 yr A/*/Mq 



where C is a parameter, and the subscript 0 represents the 
quantities when the contraction beginf] Figure 3]shows the 
fitting curve using Eq. (12 I together with our numerical re¬ 
sult. 


The protostar expands again when the next burst ac¬ 
cretion occurs and the ionizing photon emissivity decreases 
accordingly. This cycle is repeated for the following two cy¬ 
cles of quiescent and burst phases for 2000 < t (yr) < 5000. 
As the star becomes more massive, the features of the con¬ 
traction gradually diminish. The star contracts very little 
during quiescent phases when t > 5000 yr, i.e., after the 


1 Derivation of Eq. EH is as follows. Assuming that the stellar 
radius decreases roughly exponentially with a timescale given by 
Eq. |t|, i.e. the KH time, the time derivative of the stellar radius 
is 


d R* R * L*Rl 

d t t K H GMl v ' 

Using the fact that stellar luminosity is well approximated by the 
Eddington value oc M*, Eq. EH becomes 


d R* 
d t 


—const, x 


RZ 

M, 


(14) 


Because the stellar mass is almost constant during a quiescent 
phase, integrating Eq. EH gives Eq. . 
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Figure 1. Time evolution of the accretion rates (top two pan¬ 
els), stellar mass (third panel), radius (fourth panel), and ionizing 
photon emissivity (bottom panel). In each panel, the red and blue 
lines represent Models A and B (see Table ^). The evolution with 
constant accretion rates of 0.1 Mq yr —1 and 10 —3 Mq yr —1 are 
also presented with the black dashed and dotted-dashed lines re¬ 
spectively. In the bottom two panels, the lines for 10 —3 Mq yr -1 
have been horizontally shifted so that the increase of stellar mass 
matches that for 0.1 Mq yr —1 , i.e., M* = 100 Mq at the time 
of 10 3 years. The radii and ionizing photon emissivity of ZAMS 
stars are also plotted with the thick gray dashed lines, which are 
also plotted assuming the same mass variation. In the bottom 
panel the horizontal short-dashed line represents the critical ion¬ 
izing photon emissivity above which UV feedback should become 
effective with an accretion rate 0.1 Mq yr -1 (see Section |3J. 


stellar mass exceeds 500 Mq. Thus, UV feedback can be¬ 
come effective with this burst accretion model, but only for 
early evolutionary times. 

The result is understood by comparing the duration of 
the quiescent phase A£ q with the surface KH time tKH,surf 
defined in Eq. (10). For the integration range in Eq. (10), 
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Time [ yr ] 

Figure 2. The evolution of timescales with the constant accre¬ 
tion rate of 0.1 Mq yr" 1 . The black solid and dashed lines repre¬ 
sent KH timescale Ikh and accretion timescale t aC c, respectively. 
The red solid and dashed lines are the local KH and accretion 
timescales 1 k H .sn if and face,surf- 


we consider a layer covering 0.01 R* ^ r ^ R* that has 
only ~ 10 — 30 % of the total mass but 99.9999% of the 
volume. As seen in Figure [4j the surface KH time mono- 
tonically increases with time (and stellar mass), except for 
some spiky features that appear when the protostar con¬ 
tracts. The mass-radius relation deviates from Eq. § for a 
brief time period. Note that Eq. (|10[) can also be written as 


f K 11 .surf 


/ / Gm/rdm GM* 


fdl 


L*R* 


oc M, 


1/2 


(15) 


The first relation comes from the hydrostatic equilibrium 
for radiation pressure, Ts ra d ~ Pvad/p ~ Gm/r. The KH 
time increases because the stellar gravitational energy in¬ 
creases with mass as oc Ad*/R* oc M*^“ whereas the stellar 
luminosity only increases as oc M*. 

Figure [ 4 ] shows that, for t < 5 x 10 3 yr, the surface KH 
time is shorter than the duration of the quiescent phases. 
In other words, the quiescent phase is long enough for the 
protostar to lose the thermal energy trapped in the sur¬ 
face layer. The protostar then contracts. Even in such cases, 
however, the stellar contraction does not immediately follow 
the decrease of the accretion rate. After the first burst, for 
instance, the accretion rate falls down M CI at t ~ 150 yr 
(Fig. |3p , but the star begins to contract only at t ~ 400 yr. 
Figure 14] shows that the surface KH time in this epoch is 
really ~ 100 yr, which explains the above delay of the con¬ 
traction. Figure [ 3 ] shows that this delay lengthens to ~ 10 3 
yr after the second and third bursts, because the surface KH 
time increases with increasing the stellar mass as shown in 
Figure Q] 

For t > 5 x 10 3 yr, the surface KH time is longer than 
the duration of the quiescent phases. The quiescent phase is 
now too short for the star to significantly contract and the 
next accretion burst begins before the stellar surface layer 
can appreciably relax. Consequently, the evolution of the 
star after this point is almost the same as in the case with 
the constant accretion rate 0.1 Mq yr -1 . 

As discussed in Section |3.1.2| it is important to use 
tKH.surf as the thermal relaxation time for the bloated pro¬ 
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Figure 3. The same as Figure^ but for Models C and D (red 
and blue lines respectively, see Table Q. In the fourth panel, the 
black solid line represents a fitting function of the KH contraction 
for Model C (see text and footnote Q. 


tostars instead of the normal KH time Ikh since the latter 
does not consider the inhomogeneous density structure in 
the stellar interior, (see Figs. [ 5 ] and [6|. 

Since the normal KH time represents the thermal relax¬ 
ation time for the whole star, it might appear puzzling that 
Ikh is shorter than the local surface KH time fKH.surf for 
only a part of the star. The normal KH time provides the 
thermal relaxation time if the following condition is satisfied: 



GM* 

R* 


(16) 


which is not true in our case. Since most of the stellar mass 
is concentrated near the center as shown in Figures [5] and [6] 
we have 



10 - 100 


GM* 

R* 


(17) 
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Figure 4. Upper panel: We compare the KH timescales with the 
duration of the quiescent phase Ai q (horizontal short dashed line) 
in Model C. The solid and dashed lines show the time evolution 
of the local KH timescale for the bloated surface layer £kh surf 
(Eq. |l0[ ) and the normal KH timescale £kh (Eq. [t|). Lower panel: 
the evolution of the stellar radius in Model C. The red circle in the 
upper panel and vertical dot-dashed line mark the point where 
At q = £kh surf- Thus, the star can appreciably contract during 
quiescent phases only for earlier times. 



Time [ yr ] 

Figure 5. The evolution of the stellar interior structure in Model 
C. The black solid lines indicate the positions of the mass coor¬ 
dinates for 100% (stellar surface), 80%, 60%, 40% and 20% of 
the total stellar mass in descending order. The yellow regions 
denote layers with convection, according to mixing length the¬ 
ory. The brown regions show deuterium-burning layers, where the 
timescale over which deuterium mass fraction decreases by nu¬ 
clear fusion is shorter than the main-sequence lifetime. The green 
region represents a hydrogen-burning convective core, which is 
identified measuring the depletion time for hydrogen. The pink 
regions are deuterium-burning convective layers. The stellar mass 
is roughly equal to Mt for t > 10 3 yr, where M is the mean ac¬ 
cretion rate and t is time. 



Mass fraction: m/M, 


Figure 6. The mass distributions in the stellar interior in Model 
C. The three snapshots for M* = 10 2 Mg (t = 10 2 yr), 10 3 Mg 
(10 4 yr) and 4 X 10 3 Mg (4 X 10 4 yr) are presented. The latter 
two snapshots are for the largely bloated protostars. Each profile 
is plotted against the normalized mass coordinate m/M*, where 
m is the normal mass coordinate and M* is the stellar mass. The 
vertical axis also represents the normalized radial distance from 
the center, r/it*. 


depending on the stellar mass and accretion history. As a 
result, the effective thermal relaxation time for the whole 
star is now fKH.eff ~ 10 — 100 x Ikh , which is always larger 
than the surface KH time tKH.surf- 

We confirmed that tKH.eff becomes comparable to Af q 
at t ~ 1,000 yr, which does not match the epoch when 
the star stops contracting. For the bloated protostar, only 
entropy near the stellar surface determines whether the star 
contracts or not. Therefore it is reasonable to use the local 
KH time tKH.surf for the surface layer only. 


3.2.2 Variations with different burst accretion 

We can explain the results of Models A, B, and D based on 
our findings for Model C. In model B, for example, the star 
does not significantly contract during quiescent phases for 
t > 2,000 yr (Fig. [l]). Figure [7] shows the evolution of the 
KH time for Model B. We see that tKH.surf becomes equal 
to the duration of the quiescent phase At q at t ~ 2,000 yr, 
which agrees with the numerical results. 

Figures [l] and [3] show that, as expected, the star stops 
contracting earlier for quiescent phases of shorter duration. 
Unless the star contracts for a long time, the resulting in¬ 
crease of the stellar ionizing photon emissivity is not large. 
We conclude that the episodic accretion causes the stellar 
contraction and potential UV feedback only if the quiescent 
phase lasts for > 10 3 yr. 


4 DISCUSSION 

4.1 UV feedback from supermassive stars 
evolving with burst accretion 

We first discuss briefly whether UV photons from a SMS 
growing in mass due to burst accretion can ionize the sur¬ 
rounding gas to initiate strong radiative feedback. At an 
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Figure 7. Same as Figure [4] but for Model B. 


accretion rate M = 0.1 Mq yr - , the total number of hydro¬ 
gen and helium atoms infalling per second is 3 x 10 48 sec -1 . 
Assuming that the feedback becomes effective if all the 
atoms are ionized, we derive an approximate lower limit of 
ionizing photon emissivity which induces strong feedback: 
Smin = 3 x 10 48 sec -1 . This critical value is shown by the 
dashed horizontal lines in the bottom frames of Figures [T] 
and [3] In Models A and B, for which the duration of the 
quiescent phase At q is short, the ionizing photon emissivity 
is always lower than Smin (Figure [l]). By contrast, in Mod¬ 
els C and D with the longer quiescent phase duration At q , 
the stellar ionizing photon emissivity exceeds S m in when the 
star contracts in the quiescent phases. Radiative UV feed¬ 
back becomes effective in quiescent phases when At q > 10' 1 
yr, even if the mean accretion rate is 0.1 Mq yr - . 

The above estimate provides only a necessary condi¬ 
tion for appreciable UV feedback, because UV photons are 
also consumed by re-ionizing the recombining gas within an 
Hu region. Moreover, even if an Hu region appears around 
the star, its expansion is interrupted by the burst accre¬ 
tion. Suppose that the mass supply from the envelope to 
the disk continues after an Hll region is formed. This is 
likely because the Hll region first extends only in polar di¬ 
rections with respect to the disk in an early evolutionary 
stage (e.g., Hosokawa et al. J2011 ). With its mass increasing, 
the disk becomes gravitationally unstable, which then can 
cause another event of burst accretion. The stellar contrac¬ 
tion is halted and the ionizing photon emissivity suddenly 
decreases as the star bloats again (see Figs, [l] and [3|. As a 
result, the Hll region quickly disappears. Some of the gas ex¬ 
pelled by the expanding Hll region might fall back onto the 
disk before the Hll region reappears during the next stellar 
contraction. Thus, it is still uncertain to what extent the star 
can continue to grow within the context of burst accretion 
and intermittent UV feedback. Further studies are clearly 
needed to examine the overall impact of burst accretion on 
disk accretion and Hll region formation. 


4.2 Accretion histories in atomic-cooling halos 


In this paper, we have modeled the time-dependent accretion 
histories using a simple functional form with several free pa¬ 
rameters (Table 1). Despite progresses in 3D numerical sim¬ 
ulations, there remains significant uncertainty in the long¬ 
term accretion history of a supermassive star. [Latif et al.| 
(2013) follow the long-term evolution in the protostellar ac¬ 


cretion stage using the so-called sink cell technique. Their 
obtained accretion history shows some time-variability, but 
is overall rather smoother than considered in our paper. 
|Regan, Johansson fe Haehnelt| ( |2014| ) and |Becerra et al.| 
(20151 use simulations with a much higher spatial resolu¬ 


tion and find more signatures of the disk fragmentation. 
Although their simulations only follow the initial 10 — 100 
yr in the accretion stage, the results suggest that a highly 
time-dependent mass accretion can be realized in the direct 
collapse model. 


Vorobyov, DeSouza & Basu (2013) perform high reso¬ 


lution 2D simulations that follow the long-term evolution 
of a disk for normal Pop III cases. Episodic accretion oc¬ 
curs in these simulations with the mean accretion rate of 
~ 10 -3 Mq yr -1 . The mean value represents the rate onto a 
disk from the surrounding envelope. When gas accretion oc¬ 
curs through a self-gravitating disk, the accretion rate onto 
a star becomes strongly variable. When fragmentation oc¬ 
curs in the disk, clumps migrate inward to the star, and the 
accretion rate rises up to ~ 10 -2 — 0.1 Mq yr -1 . In the in¬ 
tervals between such burst events, the accretion rate falls to 
~ 10 -5 - 10 -4 Mq yr -1 . 

We also expect burst accretion with large variations 
around the mean value of ~ 0.1 Mq yr -1 in the direct col¬ 
lapse model. If the accretion rates during the quiescent phase 
falls to values M*, q < 10 -2 Mq yr -1 , the star can contract 
as suggested by our results. 

Our calculations show that the duration of the quiescent 
phase is a key parameter that determines if the star con¬ 
tracts or not. The quiescent phase will be controlled by the 
following two timescales, the effective fragmentation time 
ffrag (he. taking into account that multiple clumps can be 
produced during the fragmentation process) and the migra¬ 
tion time t mig- The former is the average time for a fragment 
to be born in a gravitationally unstable disk or l/(rate of 
clump formation). The latter is the timescale over which the 
newly formed fragments migrate inward to reach the star. 
We expect A t q ~ tf rag + t m i g . If tf rag < fmig, the duration 
of the quiescent phase will be mainly determined by t m i g . 
Conversely, if tf rag > tmi g , the duration is limited by tf rag . 

Previous numerical and analytical studies give reason¬ 
able estimates for the timescales. For example, |Vorobyov(1 


Zakhozhay & Dunham (20131 argue that, based on their 


simulations of present-day star formation, the fragmenta¬ 
tion time should be determined by the timescale over which 
the disk gets mass supply from the accretion envelope, 


f f r ;i g — 


Md 

Me 


(18) 


where Md is the disk mass and Md is the infall rate onto the 
disk. 


Inayoshi & Haiman (2014) develop an analytic model of 


the structure of an accretion disk around a SMS, and esti¬ 
mate where gravitational fragmentation occurs in the disk. 
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According to their model, there is a maximum radius R{, 
within which the disk is unstable to fragmentation. The mi¬ 
gration time for a fragment formed at Rf is then estimated 
to be 


tu 


~ 4 x 1CT yr 


(19) 


for the mean accretion rate 0.1 Mq yr -1 and the central star 
with M* < 10 4 Mq. This can be considered as the maximum 
timescale because a fragment formed at R < R{ will have 
the shorter migration time, t m j g < t m ig,m a x- 

We summarize the evolution of Af q expected from the 
above discussion as follows. In an early evolutionary phase 
for tfrag < tmig, the quiescent period Af q is limited by the 
migration time t m i g . In this case, the quiescent phase lasts 
for 10 3 yr, which is the critical value for the formation of 


an Hu region (Sec. 3.2.21, or even shorter. With the ac¬ 
cretion rate 0.1 Mq yr -1 from the envelope onto the disk, 
however, the fragmentation time tf rag exceeds the migration 
time when the disk mass reaches ~ 400 Mq. After this point, 
the quiescent phase may become longer with increasing the 
stellar mass, if the disk mass is comparable to the stellar 
mass. Unlike in the case for constant duration Af q studied 
in the present paper, the star with high mass will further 
contract during prolonged quiescent phases. The stellar ion¬ 
izing photon emissivity should be enhanced by this effect. 

Ultimately, in order to verify the above expectation, 
we need to derive realistic accretion histories realized in an 
atomic-cooling halo. To this end, multi-dimensional hydro- 
dynamic simulations would be necessary that follow the dy¬ 
namic accretion process onto growing SMSs. Future studies 
combining radiation hydrodynamics simulations and stel¬ 
lar evolution calculations should yield more realistic models 


(e.g., Hosokawa et al. 2011 Smith et al. 2012 Machida & 


Hosokawa|2013 Kuiper fc Yorke|2013 Hirano et al.|20f4l. 


5 CONCLUSIONS 


We have studied the evolution of supergiant protostars grow¬ 
ing under episodic accretion, which is expected in a self- 
gravitating circumstellar disk forming in an atomic-cooling 


halo (e.g., Vorobyov, DeSouza &; Basu|2013 

Regan, Johans- 

son & Haehnelt 2014 Inayoshi & Haiman 

2014 1 . We have 


followed the protostellar evolution with various assumed ac¬ 
cretion histories of repeating short bursts and long quiescent 
phases. In order to examine potential impact of the burst ac¬ 
cretion model, the adopted accretion histories have the same 
mean accretion rate 0.1 Mq yr -1 but different variabilities 
(Table 1). Our findings are summarized as follows: 

Burst accretion can qualitatively change the evolution 
of an accreting SMS. Although our previous calculations 
show that the stellar radius monotonically increases with 
increasing the mass under the constant accretion (Eq[9|, 
this is not always the case within the framework of burst 
accretion; i.e. the star can contract during the quiescent 
phases between the bursts. During the contraction, the stel¬ 
lar ionizing photon emissivity greatly increases because the 
effective temperature rapidly rises. For a quiescent phase of 
length At q > 10 3 yr, the ionizing photon emissivity exceeds 
3 x 10 48 sec -1 , which allows the formation of an Hu region. 
As a result, UV feedback might hinder the mass accretion 
onto the star. 


With a fixed duration of the quiescent phase Af q , stel¬ 
lar contraction can occur in an early epoch; for instance, for 
M* < 10 3 Mq with At q ~ 10 3 yr. With a longer At q , how¬ 
ever, stellar contraction will continue until the star becomes 
much more massive. This is well understood by comparing 
the duration of the quiescent phase and the surface thermal 
relaxation (or KH) time tKH.surf defined by Eq. (101. As the 
surface KH time increases with increasing the stellar mass, 
the timescale balance changes from At q > tKH.surf to the 
opposite at some point. Before this, the star contracts sig¬ 
nificantly during the quiescent phases because the bloated 
surface layer loses most of its entropy before the next burst 
begins. After the timescale inversion into At q < tKH.surf, 
the star continues to expand following roughly the mass- 
radius relation Eq. <§■ Here, the same comparison using the 
global KH time Ikh instead of tKH.surf largely overestimates 
the critical stellar mass above which the stellar contraction 
ceases. This is because the global KH time does not consider 
the inhomogeneous density structure in the stellar interior. 
Since the mass distribution is very centrally condensed when 
the star is in the supergiant stage (Figs [5]and[6|, the global 
KH time provides a poor estimate of the thermal relaxation 
time. 

Our calculations show that stellar radiative feedback 
may become important, if the quiescent phase is longer than 
~ 10 3 yr. Such long intervals between bursts of accretion are 
expected for an accretion disk around a growing SMS (pn- 
ayoshi fe Haiman|2014). If the stellar mass growth is halted 


by UV feedback, the final mass of the SMS is reduced and 
the mass of the remnant BH left behind after the star’s 
death is also reduced. Although we have adopted a simple 
functional form with parameters (Table 1) to model the ac¬ 
cretion histories, hydrodynamical simulations will provide 
more realistic estimates of the final mass. The stellar evo¬ 
lution and strength of the resulting UV radiative feedback 
will be investigated in future studies. 
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